In dit portfolio staan verschillende onderwerpen die over reproduceerbaar werken en goede workflows gaan. Voor elk onderwerp is een aparte pagina gemaakt. Op deze pagina’s staat beschreven wat ik geleerd heb tijdens mijn minor data science.
Het interperteren en verwerken van data van iemand anders is belangrijk in de data science. Er wordt hier data gebruikt van een onderzoek waar C. elegans nematodes blootgesteld aan verschillende soorten stoffen.
Belangrijkste in een data analyse is bergrijpen waar je data over gaat. Door de metadata te lezen krijg je een goed idee van waar je dataset overgaat en wat elke kopje betekent. Voor dit onderzoek willen we het effect weten van verschillende stoffen op de nematode maar ook de het effect van verschillende concentraties. Om dit te onderzoeken worden de kolommen RawData, compConcentration en expType gebruikt. In RawData staan het aantal nematodes op de plaat, in de compConcentration de concentratie van de stof en de expType de soort stof.
Om globaal een idee te krijgen van het effect worden deze variabel tegen elkaar uitget in de grafiek hieronder.
library(DT)
library(kableExtra)
library(readxl)
library(tidyverse)
tabel <- read_excel(here::here("data/CE.LIQ.FLOW.062_Tidydata.xlsx"))
tabel$compConcentration <- as.double(tabel$compConcentration)
tabel$compName <- as.factor(tabel$compName)
#x-as concentreatie in -10log
tabel$compConcentration <- ifelse(tabel$compConcentration==0,0,-log10(tabel$compConcentration))
tabel %>% ggplot(aes(x=compConcentration,y=RawData,color=compName,shape=expType))+
geom_jitter(width=0.1)+
labs(x="Concentratie (-log10) ",
y="Number of offspring")
In het experiment is de positieve controle van dit experiment is ethanol . De negatieve controle van dit experiment is de S-medium
Bij het verwerken van data is het belangrijk dat de data goed ingelezen wordt. Eén ding waar we tegen lopen bij deze dataset is dat de compConcentration ingelezen wordt als character type i.p.v. een getal. Hierdoor klopt de de x-as niet met wat je wilt zien, omdat er groepen en geen schaal weergegeven wordt. Dit wordt opgelost door het typ evan compConcentration te veranderen van een character naar een double.
De bovenstaande grafiek geeft een globaal overzicht, maar er kunnen geen conclusies uitgetrokken worden. Hiervoor is een verdere data analyse nodig. Zo’n data analyse kan er als volgt uitzien:
0,05 nulhypothese wordt aangenomen, de data is normaal verdeeld. <0,05 nul-hypothese wordt verworpen
#Genormalizeerd en gemiddelde naar 1
neg_controle <- tabel %>% filter(expType=="controlNegative")
tabel$RawData <- tabel$RawData / mean(neg_controle$RawData)
tabel %>% ggplot(aes(x=compConcentration,y=RawData,color=compName,shape=expType))+
geom_jitter(width=0.1)+
labs(x="Concentratie (-log10) ",
y="Number of offspring")
Bij de grafiek hierboven is de dataset genormaliseerd door de RawData te delen daar het gemiddelde van de negatieve controle. Dit is gedaan omdat in de concentratie en aantal nematode een grote range zit. Op deze manier kunnen hogere getallen de uitslag van analyses onterecht beïnvloeden. Door het op deze manier te normaliseren wordt de range kleiner en is er minder bias.
In dit gedeelte wordt een artikel beoordeeld op de reproduceerbaarheid. Deze wordt beoordeeld aan de hand van verschillende criteria omschreven in de onderstaande tabel. Voor de beoordeling is een artikel gekozen primair onderzoek dat beschikbaar is op PMC.
Gebruikte artikel:
Amawi KF, Alkhatib AJ. Urtica Pilulifera in Treating Pre-diabetic Rat
Model to Control the Blood Glucose, Lipids and Oxidative Stress. Med
Arch. 2020;74(3):168-171. doi:10.5455/medarh.2020.74.168-171
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7405998/
Het doel van het onderzoek is kijken of Urtica pilulifera (een plant) effect heeft op pre-diabetische ratten en ook de antioxiderende werking onderzoeken
De ratten werden ingedeeld in drie groepen van 10; een controle groep, een pre diabetische groep, en een groep met de behandeling van U. pilulifera. De ratten pre-diabetisch gemaakt kregen een hoog sucrose dieet, de controle groep een normaal dieet en de behandelde groep kreeg hetzelfde dieet als de pre-diabetische ratten met U. pilulifera extract geïnjecteerd. Na 30 dagen werden bloed samples afgenomen en getest op glucose, triglycerides, cholesterol, GSH, TAC en MDA
Uit het onderoek bleek dat glucose, triglyceride en MDA niveaus in de pre-diabetic groep significant verhoogd waren en significant verlaagd in de U. pilulifera groep. GSH en TAC was significant hoger in de U. pilulifera ten opzichte van de pre-diabetische groep. Er zat geen significant verschil in cholesterol niveau in de groepen.
In de tabel staat of het artikel aan de verschillende criterium voldoet.
| Criteria | Answer |
|---|---|
| Study purpose | Yes |
| Data availability | No |
| Data location | Data location not stated |
| Study location | Yes |
| Author review | Author listed but did not fill out contact |
| Ethics statement | No |
| Funding statement | Yes |
| Code availability | No |
Het blijkt dat het artikel maar aan drie van acht criteria voldoet. Hoewel het artikel goed te lezen was is deze dus toch niet goed reproduceerbaar.
Net is er gekeken reporduceerbaarheid van een primair onderzoek. In dit onderdeel wordt er gekeken naar de reporduceerbaarheid van (R) code.
Data en code van onderstaand onderzoek gebruikt: Strozza C, Myrskylä M. Monitoring trends and differences in COVID-19 case-fatality rates using decomposition methods: Contributions of age structure and age-specific fatality. PLoS One. 2020 Sep 10;15(9):e0238904. doi: 10.1371/journal.pone.0238904. PMID: 32913365; PMCID: PMC7482960. https://pubmed.ncbi.nlm.nih.gov/32913365/
Link naar de code: https://osf.io/g7vjd/
Het script zelf was goed te lezen. Er waren 4 scripts die elkaar opvolgde, door de getallen in de bestandsnaam was het makkelijk te zien welke als eerst uitgevoerd moest worden. In de scripts zelf stonden genoeg comments om te begrijpen waar elk stuk code voor diende en ziet de code zelf er netjes uit. Ik geef de leesbaarheid dus ook een 5/5
De code was ook goed reproduceerbaar. Er was alleen één probleem wat handmatig opgelost moest worden, dit probleem en hoe deze opgelost is staat beschreven hieronder. Hierdoor krijgt de reproduceerbaarheid een 4/5, want nadat dit probleem opgelost was, konden alle tabellen makkelijk gemaakt worden.
### Monitoring trends and differences in COVID-19 case fatality ##############
### rates using decomposition methods: A demographic perspective ##############
### Last updated: 2020-07-16 09:27:20 CEST
### Contact:
### riffe@demogr.mpg.de
### acosta@demogr.mpg.de
### dudel@demogr.mpg.de
### Get data ##################################################################
# Required packages
source(("R/00_functions.R"))
# URL + filename
url <- 'https://osf.io/wu5ve//?action=download'
filename <- 'Data/Output_10.csv'
# Load data
GET(url, write_disk(filename, overwrite = TRUE))
dat <- read_csv(filename,skip=3)
Bij het uitvoeren van de code hierboven kwam de onderstaande foutmelding:
In de foutmelding
is een html bestand te zien. Dit hoort natuurlijk niet bij een
read_csv() functie. In de code wordt het bestand van een online
directory gedownload. De URL staat in de code en wanneer deze bezocht
werd stond er dit bericht:
Error op de webpagina
Het bestand bevindt zich niet meer op de plek waar de URL naar verwees en kon dus ook niet gedownload worden waardoor de html van de pagina werd overgenomen. Er was echter wel een menu aan de zijkant aanwezig waar het mogelijk is om bij de github van het project te komen. Hier vond ik het output bestand wat nodig is voor de code en heb deze handmatig gedownload en toegevoegd aan de /data directory. Nadat deze in de directory stond kon het script en de rest van de scripts vlekkeloos uitgevoerd worden. Hieronder staat het script met de aanpassing. De rest van het script kan gevonden worden op de bovenstaande link. De output van het script zijn 6 verschillende excel bestanden met data over het verloop van COVID-19 in 6 verschillende landen.
### Monitoring trends and differences in COVID-19 case fatality ##############
### rates using decomposition methods: A demographic perspective ##############
### Last updated: 2020-07-16 09:27:20 CEST
### Contact:
### riffe@demogr.mpg.de
### acosta@demogr.mpg.de
### dudel@demogr.mpg.de
### Get data ##################################################################
# Required packages
source(("R/00_functions.R"))
# URL + filename
filename <- 'data/Output_10.csv'
# Load data
dat <- read_csv(filename,skip=3)
### Edit data (select countries, etc.) ########################################
# Lists of countries and regions
countrylist <- c("China","Germany","Italy","South Korea","Spain","USA")
region <- c("All","NYC")
# Restrict
dat <- dat %>% filter(Country %in% countrylist & Region %in% region)
# Remove Tests variable
dat <- dat %>% mutate(Tests=NULL)
# Drop if no cases/Deaths
dat <- na.omit(dat)
### Save ######################################################################
write_csv(dat,file="Data/inputdata.csv")
Het organizeren van data is belangrijk voor reproduceerbaarheid en algemene goede workflow. Hieronder staat een voorbeeld van een map met verschillende projecten. Elke map heeft dezelfde indeling, bij de data en scripts zijn README bestanden toegevoegd waarin metadata staat over het project en de bestanden. Ik vind het zelf fijn dat wanneer er veel van één soort bestand is (zoals bij bam en fastq bestanden) dat deze in een aparte map staan. De rest bestanden staan gewoon in de /data.
## data/daur2
## +-- metagenomics_formatief
## | +-- data
## | | +-- fastq
## | | +-- HU_waternet_MOCK2_composition.csv
## | | \-- README.md
## | +-- metagenomics_formatief.html
## | +-- metagenomics_formatief.Rmd
## | +-- metagenomics_reader.Rproj
## | \-- README.md
## +-- metagenomics_reader
## | +-- data
## | | +-- fastq
## | | +-- HU_waternet_MOCK1_composition.csv
## | | \-- README.md
## | +-- metagenomics_reader.html
## | +-- metagenomics_reader.Rmd
## | +-- metagenomics_reader.Rproj
## | \-- README.md
## +-- rna_seq_airway
## | +-- data
## | | +-- bam
## | | +-- counts
## | | +-- fastqc
## | | +-- fastqc_output
## | | +-- ipsc_sampledata.csv
## | | \-- README.md
## | +-- README.md
## | +-- rnaseq_airway.html
## | +-- rnaseq_airway.Rmd
## | \-- rnaseq_airway.Rproj
## +-- rna_seq_ipsc
## | +-- data
## | | +-- bam
## | | +-- counts
## | | +-- fastqc
## | | +-- fastqc_output
## | | +-- ipsc_sampledata.csv
## | | \-- README.md
## | +-- README.md
## | +-- rnaseq_ipsc.html
## | +-- rnaseq_ipsc.Rmd
## | \-- rnaseq_ipsc.Rproj
## \-- rna_seq_onecut
## +-- data
## | +-- bam
## | +-- counts
## | +-- fastqc
## | +-- fastqc_output
## | +-- onecut_sampledata_OC2.csv
## | \-- README.md
## +-- README.md
## +-- rnaseq_onecut.html
## +-- rnaseq_onecut.Rmd
## \-- rnaseq_onecut.Rproj
Hieronder staat een voorbeeld van de readme die ik in de repo van dit portfolioo heb staan. Ik beschrijf kort wat er in de datasets staat en waarvoor ze gebruikt. Ook staat er waar ik ze vandaan gehaald heb.
README voorbeeld bij data bestanden
Foto van mij
0612345678
lisa.kuipers26@hotmail.com
HEMA Rhenen | 2015 – heden
Kassa medewerker, horeca medewerker
Rembrandt college Veenendaal | 2012 - 2017
Hogeschool Utrecht | 2019 – Heden
HBO Biologie en medisch laboratoriumonderzoek (life science)
Nederlands
Engels
Microsoft office
R
Bash
Python
SQL
High trougputscreening is een techniek die gebruik wordt naar het onderzoeken van effecte van bepaalde medicijnen op kankercellen. Met deze technieken worden veel verschillende medicatie getest op cellen waaruit veel data uit één keer ontstaat(Carnero 2006) De opdracht werd gegeven door het prinses maxima center. Die zich focused op het onderzoek naar kinderoncologie(“Ons verhaal” n.d.) Onder gedaan naar leukemie (Brivio et al. 2021)
Door gebruik te maken van shiny kan de programmeertaal R gebruikt worden om apps te maken (“Shiny” (n.d.))
“Survival Rate for Children with Brain Tumors Improved Since 1990” (n.d.)
library(dslabs)
library(tidyverse)
#laden van gapminder
gapminder_df <- gapminder %>% as.data.frame()
#tidy maken van data functie
tidy_func <- function(path, ziekte){
read_df <- read.csv(path, skip=11) %>% as.data.frame
col_nam <- colnames(read_df)
df_tidy <- read_df %>% pivot_longer(cols=col_nam[-1], names_to="country", values_to=paste0(ziekte,"_activity")) %>%
separate(Date, into = c("year", "month","day"),convert = TRUE)
}
dengue_tidy <- tidy_func("data/dengue_data.csv","dengue")
flu_tidy <- tidy_func("data/flu_data.csv","flu")
#opslaan als RDS
saveRDS(dengue_tidy, file = "data/dengue.rds")
saveRDS(flu_tidy, file = "data/flu.rds")
saveRDS(gapminder_df, file = "data/gapminder.rds")
#opslaan als CSV
write.csv(gapminder_df,"data/gapminder.csv")
write.csv(dengue_tidy,"data/dengue_tidy.csv")
write.csv(flu_tidy,"data/flu_tidy.csv")
SQL code voor het creëren van een database:
CREATE DATABASE workflowsdb;
library(DBI)
#Connectie opzetten
con <- dbConnect(RPostgres::Postgres(),
dbname = "workflowsdb",
host="localhost",
port="5432",
user="postgres",
password="Datascience")
#Aanmaken van tabellen in databse
dbWriteTable(con, "flu", flu_tidy)
dbWriteTable(con, "dengue", dengue_tidy)
dbWriteTable(con, "gapminder", gapminder_df)
#Jaren verwijderen die niet aanwezig zijn in flu en dengue
gapminder_minyear <- gapminder_df[gapminder_df$year >= min(flu_tidy$year), ]
gapminder_maxyear <- gapminder_minyear[gapminder_minyear$year <= max(flu_tidy$year), ]
#Landen weghalen die niet en flu en dengue aanwezig zijn
landen <- c(col_dengue[2:11], col_flu[2:30],recursive=TRUE)
gapminder_clean <- subset(gapminder_maxyear, country %in% landen)
#Oude tabel verwijderen en nieuwe tabel invoegen
dbRemoveTable(con, "flu")
dbWriteTable(con, "gapminder", gapminder_clean)
De ondesrtaande code is sql code wat gebruikt wordt om een grote tabel te creëren waar alle data in staat. Er worden left joins gebruikt omdat niet alle landen in beide de dengue dataset en flu dataset voorkomen. Deze informatie van bepaalde datasets kan wel interessant zijn voor het analyseren van data en wordt met de left join een NULL ingevoerd bij de missende data. Voor de rest worden alle kollomen van gapminder gekozen, de dag en maand van flu en de activity van beide flu en dengue. ook worden er joins op year, month en day gedaan. Dit is zodat ze gematched worden in de tabel en laad tijd verminderd bij het sorteren. Bij het sorteren wordt eerst op het jaartal gesorteerd en daarna op het land om het zo overzichtelijk te houden.
SELECT
gapminder.*,
flu."day",
flu."month",
flu.flu_activity,
dengue.dengue_activity
FROM
gapminder
LEFT JOIN flu
ON gapminder."year" = flu."year" and gapminder.country = flu.country
LEFT join dengue
ON dengue."year" = flu."year" and dengue."month" = flu."month" and dengue."day" = flu."day" and dengue.country = flu.country
ORDER BY gapminder."year", gapminder.country ;
Deze query wordt uitgevoerd in het programma beaver, hieronder een foto van het dashboard
FOTO !
Door de volgende tekst toe te voegen boven bovenstaande query wordt er een tabel gemaakt genaamd all_data
CREATE TABLE all_data
AS
(...)
Op deze afbeelding is te zien dat de tabel toegevoegd is aan de database
FOTO !
#Data uit DB ophalen
all_data <- dbReadTable(con, "all_data")
#Dataset opslaan in data
saveRDS(all_data, file="data/all_data.rds")
#Data laden voor analyse
all_data <- readRDS("data/all_data.rds")
#Functie voor life exptancy tegen andere kolom
lifeexp_func <- function(dataset, kolom){
means <- dataset %>%
group_by(year) %>%
summarise(naam=mean(!! sym(kolom),na.rm=TRUE),
"mean_life_expectancy"=mean(life_expectancy,na.rm=TRUE))
colnames(means) <- c("year",paste0("mean_",kolom),"mean_life_expectancy")
means
}
#Flu en dengue tegen life expactancy
mean_flu_year <- lifeexp_func(all_data,"flu_activity")
mean_dengue_year <- lifeexp_func(all_data,"dengue_activity")
mean_flu_year %>% ggplot(aes(x=year,y=mean_life_expectancy))+
geom_line()
mean_flu_year %>% ggplot(aes(x=year,y=mean_flu_activity))+
geom_line()
mean_dengue_year %>% ggplot(aes(x=year,y=mean_dengue_activity))+
geom_line()
#country_data <- all_data %>% filter(year==2008) %>%
# group_by(country) %>%
# summarise(mean(flu_activity), max(population),max(gdp)) %>% na.omit()
#colnames(country_data) <- c("country","flu_activity","population","GDP")
#country_data$population <- factor(country_data$population, levels=c(arrange(desc(country_data$population))))
library(RColorBrewer)
#Summarise van data per year
flu_dengue <- all_data %>% group_by(month, year) %>%
na.omit() %>% filter(year!=2002) %>%
summarise(mean(dengue_activity, na.rm=TRUE),mean(flu_activity, na.rm=TRUE))
colnames(flu_dengue) <- c("month","year","dengue_activity","flu_activity")
#Functie voor plot van verloop activity
activity_years_func <- function(dataset,activity){
dataset %>% ggplot(aes(x=month,y=dataset[[activity]], group=year, color=factor(year)))+
geom_line()+
scale_x_continuous(name="Month", breaks=(1:12), limits=c(1, 12))+
scale_color_brewer(palette = "Set1")+
theme_minimal()
}
activity_years_func(flu_dengue,"flu_activity")
activity_years_func(flu_dengue,"dengue_activity")
Parameters in yaml headers kunnen gebruikt worden om je code dynamischer te maken. Het si een snelle manier van variabele veranderen van bijvoorbeel dhet filteren van een dataset.
Hieronder staat een voorbeeld van hoe de yaml header eruit ziet met verschillende parameters. Voor de analyse verderop hebben we dezelfde paramteres gebruitk als in de afbeelding.
Afbeelding 1: Parameters in de yaml header.
Afbeelding 2: Parameters bij knitten.
Hieronder staat de code van het filteren van de de EDCC_daily dataset found at https://www.ecdc.europa.eu/en/covid-19/data. In de datset staan het aantal gevallen en sterfgevallen van COVID-19 uit de periode van 2020 t/m 2022. De parameters die boven staan zijn gebruikt in het filteren. Wanneer iemand van 1 jaar maand of dag de data wilt weergeven kan bij de end en start parameters dezelfde datum ingevoerd worden. Ikzelf heb ervoor gekozen om de hele periode te nemen.
library(tidyverse)
library(ggplot2)
library(params)
library(stringr)
data_table <- read_csv("data/EDCC_daily.csv")
data_filtered <- data_table %>% filter(year >= params$start_year & year <= params$end_year & countriesAndTerritories == params$country & day %in% c(params$start_day:params$end_day) & month %in% c(params$start_month:params$end_month))
data_filtered$date_range <- paste0(data_filtered$year,sprintf("%02d",data_filtered$month),sprintf("%02d",data_filtered$day))
plot_func <- function(dataset,condition){
dataset %>% ggplot(aes(x=date_range, y=dataset[[condition]]))+
geom_point()+
scale_x_discrete(breaks=data_filtered$date_range[seq(1,length(data_filtered$date_range),by=50)])+
theme(axis.text.x = element_text(angle = 90))+
labs(x="Date",
y=condition,
title = paste("Amount of COVID-19",condition,"from",data_filtered$dateRep[length(data_filtered$dateRep)],"to",data_filtered$dateRep[1],"in Belgium"))
}
plot_func(data_filtered,"cases")
Figuur 1: Aantal COVID-19 cases van 01/03/2020 t/m 23/05/2022 in belgië. Aantal cases staat op de y-as en de datum als yyyy-mm-dd op de x-as
plot_func(data_filtered,"deaths")
Figuur 2: Aantal COVID-19 deaths van 01/03/2020 t/m 23/05/2022 in belgië. Aantal cases staat op de y-as en de datum als yyyy-mm-dd op de x-as
Wat interessant is om te zien bij deze figuren dat eigenlijk het aantal doden aan het begin het hoogste lag terwijl het aantal cases niet heel hoog waren. Bij de eerste piek van cases is bij het aantal dode ook een piek te zien ronde dezelfde datum. Daarna zwakt het aantal doden vergeleken met de cases af. Of dit komt doordat mensen met een verzwakt imuunsysteem sneller zijn overleden of omdat vaccinaties werken is natuurlijk niet direct op te maken uit deze data, maar het is wel een interessant verband.